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In this work we study the spreading dynamics of tiny liquid 
droplets on solid surfaces in the case where the ends of the 
molecules feel different interactions with respect to the sur- 
face. We consider a simple model of dimers and short chain- 
like molecules that cannot form chemical bonds with the sur- 
face. We use constant temperature Molecular Dynamics tech- 
niques to examine in detail the microscopic structure of the 
time dependent precursor film. We find that in some cases 
it can exhibit a high degree of local order that can persist 
even for flexible chains. Our model also reproduces the ex- 
perimentally observed early and late-time spreading regimes 
where the radius of the film grows oc t 0,5 . The ratios of the 
associated transport coefficients are in good overall agreement 
with experiments. Our density profiles are also in good agree- 
ment with measurements on the spreading of molecules on 
hydrophobic surfaces. 

68.10. Gw, 61.20.Ja, 05.70.Ln 



I. INTRODUCTION 

Since the pioneering spreading experiments of micro- 
scopic liquid droplets on surfaces by Heslot et al. |Q, 
there has been increasing interest in phenomena occur- 
ing at microscopic lengthscales during spreading. The 
experiments of Refs. reveal that both the molec- 

ular structure of the liquid and the type of substrate 
used influence density profiles of the droplets. For exam- 
ple, thickness profiles of tetrakis (2-ethylexoxy)-silane 
and polydimethylsiloxanc (PDMS) droplets on a silicon 
wafer exhibit strikingly different shapes under spread- 
ing [Q. Tetrakis exhibits clearly observable dynamical 
layering, while the spreading of PDMS proceeds by a 
quickly evolving precursor layer of one molecular thick- 
ness. Furthermore, the experiments of Valignat et al. |j| 
and Cazabat et al. || address the important role of sur- 
face grafting, and asymmetrical surface interactions in 
particular, on the spreading dynamics. For example, in 
Rcf. U the density profiles for trisiloxane polyoxyethy- 
lene droplets spreading on a hydrophilic surface resemble 
a "sand pile" , whereas the same liquid forms a very com- 
pact bilayer on a hydrophobic surface. 

Despite drastic differences in the density profiles, an 
interesting feature in the experiments of Refs. is 



that the time dependence of the radius of the precursor 
film r(t) follows "diffusive" behavior r(t) ~ i 5 for all 
times measured. Moreover, the experiments of Valignat 
et al. [jsj report two distinct "diffusive" regimes compris- 
ing a rapid early-time region followed by a considerably 
slower late-time one. Typical estimates for the ratio be- 
tween the corresponding early-time and late-time trans- 
port coefficients are of the order of 100. The emergence 
of the early-time regime can most simply be explained 
by assuming that the flux onto the surface is constant, 
i.e. dN p (t)/dt = const., where N p (t) is the number of 
molecules in the effectively 2D precursor film H. From 
this it follows that r(t) ~ y/N p (t) ~ t - 5 . The late-time 
behavior is due to crossover towards 2D diffusion that 
takes over in the submonolayer regime QJ^]. 

Motivated by these experimental discoveries a num- 
ber of theoretical models have been proposed. However, 
progress has been rather moderate. Analytic theories to 
date deal with dynamical layering only ||Jl^] . Abraham 
et al. || considered a solid-on-solid (SOS) model of a 
layered droplet in rectangular geometry. The width of 
the precursor film was found to evolve linearily in time, 
implying that the effective flux of the liquid into the pre- 
cursor film is indeed constant in time. De Gennes and 
Cazabat [|l0| considered a model in which the layers have 
already formed. This model gives simple relations for the 
time dependence of the radii of the layers. In particular, 
for a completely layered droplet, the precursor film de- 
velops approximately "diffusively" in time JTT[ . However, 
both models treat the liquid as structureless and there- 
fore cannot be applied to studying the effects of molec- 
ular structure or details of interactions on the spreading 
dynamics. 

For both coarse-grained and more microscopic models, 
a number of computer simulations have been performed 
]l2| |i~9[ to study the dynamics of spreading. In particular, 
using Molecular Dynamics (MD) simulations it was con- 
cluded in Ref. jl6j that both the chain-like nature of the 
molecules and the chain-surface interactions can signifi- 
cantly influence the structure of the precursor layer. Us- 
ing a cylindrical droplet geometry, Refs. 0,^6) reported 
an "almost linear" early time regime for the width w(t) 
of their precursor film (again indicating a constant flux), 
followed by a late-time diffusive region. Most recently, 
the t ' 5 behavior was observed in another MD simulation 
|l9| . Qualitatively, the two time regimes have been seen 
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in Monte Carlo simulations, too |T^,|T^]. 

In the present work, our aim is to employ MD simula- 
tion techniques to study a particularly interesting aspect 
of droplet spreading, namely the case where the fluid 
molecules can feel asymmetric interactions with respect 
to the surface. We will concentrate mostly on the mi- 
croscopic structure of the precursor film, its time depen- 
dence and the quantitative evaluation of the associated 
transport coefficients. Our work is motivated by recent 
experiments on such systems |^||, as well as the practi- 
cal importance of the molecular structure of thin layers. 
In particular, a common way of controlling the surface 
energetics of a solid substrate is to use grafted molecules 
that adsorb on iLsometimes forming chemical bonds and 
brushed layers [p0| . Such surfaces have important appli- 
cations in e.g. coating and lubrication. Another inter- 
esting class of systems comprises amphiphilic molecules 
such as detergents where a strong asymmetry of interac- 
tions causes layered structures to form Spreading 
dynamics of such molecules is then of particular interest 
in trying to understand how these layers form, and how 
well ordered they are. 

The outline of this paper is as follows. First, we intro- 
duce the model in detail. Then we present results for the 
spreading dynamics of molecules consisting of two, four, 
and eight units with an asymmetry of interaction with re- 
spect to the underlying surface. Our results reveal that 
in some cases the emerging precursor layer may exhibit 
a high degree of local order that persists even for the 
longer chains. This is reflected in the corresponding den- 
sity profiles, some of which bear a close resemblance to 
the experimental ones §. We also examine the time de- 
pendence of the radius of the precursor layer, and find the 
two different "diffusive" regimes in agreement with previ- 
ous studies 0,|l6| and experiments ||. We furthermore 
evaluate the corresponding transport coefficients and find 
that their ratio is in good qualitative agreement with the 
experiments. Finally, we briefly discuss the influence of 
the choice of thermostat in our model. A brief account 
for the results for dimers has been previously given in 
Ref. H). 

II. MODEL 

A. Interactions 

Our MD model is analogous to the one introduced in 
Refs. 0,^6|. The n-mer molecules consist of n Lennard- 
Jones (LJ) particles. Within the chain, the LJ interac- 
tion between them is purely repulsive, i.e. of the form 
Vl r j tra (r) — 4e/(cr//r) 12 to prevent spatial overlap. The 
potential parameters are <jf = 2.3 A for the width and 
tf = 0.1703 eV for the depth, respectively. Additionally, 
the particles are interconnected by a very rigid but ori- 
cntationally isotropic harmonic oscillator pair potential 
V c = ifc(r-r ) 2 , where k = 10000e//oj and r = 2 1 / 6 ct / 



so that the chains do not easily stretch. There is also 
an angle dependent potential Vg — eg(cosd + 1), for 
n > 2. For our studies we examine two cases, namely 
eg = 10e/ (rather stiff chains) and eg — (completely 
flexible chains). There is no torsion dependent potential 
within a chain since we consider linear chains. 

Interchain interactions are modeled by the following 
pairwise LJ interaction between n LJ monomers: 

V LJ (r) = 4e f [(^r - (^f] . (1) 

The substrate on the other hand is modeled by a flat con- 
tinuum LJ material; it is thought to be homogeneous and 
its unit volume = 1. The total substrate interaction 
is obtained by integrating the LJ potential over the half 
space z < with the result 

W = -^ + f> ( 2 ) 

where 

A = {27r/3) Ps e l <jS and B = (4 7 r/45) j0s e l( T. i 12 , (3) 

and where p s denotes the number density of particles 
comprising the substrate. Different chain-surface inter- 
action parameters employed in this study are presented 
in Table I and shown in Fig. 1. 

The asymmetrical nature of the chain-surface inter- 
actions comes about through the choice of different set 
of surface interaction parameters for the grafted end as 
compared to the other monomer units along the chain. 
The "grafted" end interaction is set to Vj in every case 
studied in this work, whereas for dimers we employ Vjj 
(ordinary case) and Vjjj (shifted case), for tetramers Vjy 
(ordinary case) and Vy (shifted case), and for octamers 
Vjy. These potentials have been constructed in such a 
way that in the ordinary case individual chains tend to 
lie parallel to the surface while in the shifted cases they lie 
perpendicular to it. It should be noted that the chains in 
our model are not allowed to form chemical bonds with 
the surface. 



B. Choice of physical units 

The physical units are determined by the Hamaker 
constant Ah of the substrate and by the number den- 
sity of substrate atoms. In our units, we have fixed 
p s = 1.0 A . The Hamaker constant Ah of the sub- 
strate fixes the effective bond length hi . This can be seen 
as follows: the Hamaker constant between two materials 
is defined by 

A H = 4tt 2 e s a 6 s p s p f , (4) 

where p/ denotes the number density of molecules in the 
fluid [pTJ . We have fixed Ah through the choice e s = 
2.8 x 10~ 20 J and cr s — 1.25 x 10~ 9 m and requiring that 
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it is realistic, i.e. A H ~ 10~ 18 J @. This fixes the 
number density of the fluid to be pf ~ 10 24 m -3 . On 
the other hand, pf ~ and hence the effective bond 
length bi ~ 10~ 8 m = 100 A. This justifies the use of a 
flat, continuum substrate in our model studies. 

For an LJ system there is a typical time scale, which 
is fixed by the choice of e/, 07, and the mass to. In 
our bare units, the mass of a monomer is mb = 63.5 
amu, e/ = 0.1703 eV and 07 = 2.3 A. The combination 
that yields a quantity with the dimension of time is r c = 

{rrib(J 2 j)/€f = ^/nrib/ef a / which in our bare units is « 

5 x 10 -13 s. In our simulation algorithm, we have chosen 
the time step to be 0.01r c . To obtain physical units we 
use the bond length bi ~ 100 A which scales 07, and 
set the physical mass of our effective monomers to be 
a realistic value of m = 10 5 amu. Using these values 
to scale t c gives the time step in physical units to be 
t TM . = 7.7 x 10" 13 s. 

C. Choice of thermostat 

The dynamics of the system with a Nose-Hoover (NH) 
thermostat is described by the usual equations of motion 



^ = _ Vl u - m , (6) 

and 

^E^-ww, (7) 

i 

where Nf is the number of degrees freedom, kT s the tem- 
perature for the thermostat, r)(t) a time-dependent fric- 
tion coefficient, and r = 2.0 x 10~ 14 s is a relaxation 
time. If we were to study the equilibrium properties 
of our model, we note that this choice of r would re- 
move any uncanonical temperature fluctuations due to a 
hidden Toda demon [E3|. The equations of motion are 
solved using modified velocity Verlet algorithm (see e.g. 
Ref. The simulations are performed at temperature 

kT = 0.8e/ which is well above the triple point of an LJ 
fluid (l(| . At the end of this work, we will also briefly 
discuss the influence of choosing a local thermostat based 
on Langevin dynamics to our results. 

D. Construction of the initial configuration 

We use the cylindrical geometry of Refs. |l4|,[l6| an d 
construct an initial ridge-shaped droplet with periodic 
boundary conditions along the direction of the ridge 



which is denoted by y. The length of the cylinder is 
« 10 6; for dimers, ~ 16 6; for tetramers and « 38 6; for 
octamers. Spreading takes place in the x direction which 
lies perpendicular to the ridge. This choice of geometry 
is for computational convenience, and translating our re- 
sults to the true 3D case is straightforward. 

Constructing a proper initial configuration is com- 
plicated by long relaxation times of the chain-like 
molecules. To overcome this, in the beginning the lo- 
cations of the end-groups of the n-mers are chosen ran- 
domly. Then the chains are formed in such a way, that 
they are bent 90 degrees at each joint, while their direc- 
tions are random. The initially rather sparse system is 
compressed to find the minimum of the internal energy. 
Then the droplet is allowed to equilibrate in the follow- 
ing way: the temperature of the system is set to 0.1 kT s , 
where kT s = 0.8 e/ denotes the actual simulation temper- 
ature. Then the temperature is raised to 1.5 kT s by con- 
tinuously adjusting the temperature of the thermostat. 
All this time the system is allowed to evolve without the 
surface interaction. After this the temperature is lowered 
to kT s in the same way; this should result in a configura- 
tion that is closer to the actual equilibrium configuration 
than the initial with. We have qualitatively confirmed 
this by estimating the corresponding free energy differ- 
ences . The time scales used in constructing the ini- 
tial configuration are comparable to the actual spreading 
simulation. After the system has been equilibrated, the 
substrate potential is "switched on" and spreading can 
take place. 

E. Quantities calculated 

One of the main advantages of the MD method is that 
the spreading dynamics can be followed in real time, and 
detailed data on the configurations are available. To this 
end, we have calculated the following quantities: 

(i) The width of the precursor film w(t) is a mea- 
sure of the average horizontal extent of the droplet 
in the x direction. In our geometry the spreading 
takes place in the x direction only and hence sim- 
ple geometrical arguments give w(t) « N p (t)/(pL), 
where N p (t) is the number of particles in the pre- 
cursor film, p denotes the average number of parti- 
cles per unit area in it and L denotes the length of 
the droplet along the y axis. To convert this to the 
spherically symmetric 3D case we assume that the 
flux of particles to the precursor film cx t a where 
awl. This immediately implies, that w(t) * At a 
for a compact layer. It follows then that for our 
geometry 

N p (t) w ApLt a . (8) 

On the other hand, the radius of the precursor film 
of a fully spherical droplet that spreads radially, is 
given by 
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r(t) « sjN p {t)/{pir) . 

Combining Eqns. (g) and (||) we obtain 

r(t) Pa y/ApLV*/{jm) 
= yfALj^ t a ' 2 



(9) 



(v) 



(10) 



where the associated early-time "diffusion" (trans- 
port) coefficient is defined by D e = LA/tt. 

(ii) The density profile is of particular interest be- 
cause it can be directly measured in the experi- 
ments. In our model it is defined as the average 
number of particles in a box of size Ax Ay taken at 
any fixed y. 

(iii) The pair correlation function characterizes the 
degree of order in the system. It is defined as the 
following configuration average p4fl: 



g(r,t) = (N/V)- 2 (J2Y, S ^ Cm W f i 



(11) 



where N denotes the number of molecules in the 
system of volume V and t f i an is the position of the 
centre-of-mass of the i th molecule. In particular, 
we have calculated g(r,t) within the effectively 2D 
precursor film, where N = N p (t) and V = A p (t), 
the latter denoting the time-dependent area of the 
precursor film. It should be noted that since the 
system is not translationally invariant in the x di- 
rection, g(r) starts deviating from its asymptotic 
values for large r, for which g(r)—tl when r— >oo. 

(iv) The distribution of orientations for neigh- 
bouring chains d Q is defined as follows: to every 
chain % a vector B4 is assigned as B4 = fi +n i — 7$. 
Here fi denotes the position of the grafted end of 
the chain i, and ri+ n i denotes the position of the 
other end of the chain. The normalized distribution 
is defined as follows: 



d (9,t) 



N o (0,t) 



(12) 



where N a (9, t) denotes the number of pairs of chains 
oriented at an angle 9 relative to each other at 
time t. Only those pairs of chains whose centres of 
masses lie within a sphere of radius 2.5 at are con- 
sidered. This quantity indicates the spatial orien- 
tational correlations of neighbouring chains. This 
method is accurate when we use rigid chains, but 
in the case of flexible ones it gives only an approxi- 
mate picture of the degree of orientational correla- 
tions between neighbouring chains. 



Distribution of bond angles db is a time- 
dependent quantity that keeps track of bond angles 
at given times. The bond angle is defined to be the 
angle between consecutive bonds in a molecule and 
it is calculated by taking the dot product of the 
associated bond vectors. Bond vector is defined as 
the vector between consecutive monomers in the 
chain. Distribution of bond angles db is defined as 
the following normalized histogram: 



d b (e,t) 



N b (6,t) 



(13) 



where Nb(9,t) denotes the number of bond angles 
with angle 9 at time t. It is calculated for the initial 
configuration and for a configuration taken at late 
stages of spreading. It is used in analyzing whether 
completely flexible chains become effectively stiffcr 
in the course of spreading or not. This tendency 
would be revealed by a distribution profile peaked 
sharply around 9 w 0°. This quantity is calculated 
for tetramers and octamers only. 



III. RESULTS 

First we consider short, rigid molecules with two 
monomer units. These are called dimers, and we will 
present results for two different sets of surface interaction 
parameters. Then we consider longer flexible chains, and 
we present results for chains with four monomer units 
(tetramers) and with eight units (octamers). 



A. Dimers 

In this section we present complete results for short 
and rigid molecules (see Ref. [^2) for a brief summary). 
We have studied two different cases, namely one in which 
the equilibrium orientation of an individual chain is par- 
allel to the surface (the ordinary case) and the other in 
which the molecules prefer to lie perpendicular to the 
surface (the shifted case). The shifted case can be con- 
sidered to be an effective model for rigid molecules with 
one end hydrophobic and the other hydrophilic relative 
to the surface The grafted end interacts with the 
surface with potential Vj while in the ordinary case the 
other end has Vjj and in the shifted case Vttt. 

The ordinary case 

Figs. 2(a) and (b) show a sequence of snapshots from 
a typical evolution of the droplet for TV = 1525 ordinary 
dimers during spreading as seen along the axis of the 
ridge. The holes in the initial configuration are due to 
density fluctuations. The initial configuration is charac- 
terised by a disordered and liquid-like structure. This 
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is revealed by the pair correlation function and the dis- 
tribution of orientations for neighboring dimers for this 
geometry. 

After switching on the surface attraction it can be seen 
that the dimers in the middle of the droplet that have 
not yet come into contact with the surface are mainly 
oriented perpendicular to the surface with the grafted 
end pointing downwards due to stronger surface attrac- 
tion. The molecules that are on the surface, on the other 
hand, behave quite differently. It can be seen that the 
precursor film is very disordered except very close to the 
edges of the droplet. This is due to the high density of 
dimers near the center of the droplet which effectively 
prevents them from attaining their equilibrium orienta- 
tions relative to the surface. At the edges of the droplet 
the dimers have enough room to lie flat on the surface. 
The final configurational stage is a thinning monolayer 
of molecules lying flat on the surface and exhibiting dif- 
fusive motion. 

To quantify these observations we have calculated the 
pair correlation function of the center of mass of the 
molecules within the precursor film. This case shown 
in Fig. 3 reveals that there is only weak short-range or- 
der characteristic of liquids. In order to characterize the 
degree of long-range orientational correlations we have 
calculated the distribution of orientations for neighbour- 
ing dimers at a late stage of spreading. It reveals that the 
overall orientational correlation between well-separated 
dimers is weak. This again is consistent with conclu- 
sions made above about the structural order within the 
droplet. 

Fig. 4 shows the evolution of the density profile of 
the droplet taken at three different time steps. It can 
be seen that the profiles are fairly smooth and rounded; 
the peaks are due to dimers that have not yet come into 
contact with the surface. At later times a step develops 
at the edge of the film where dimers tend to lie flat on 
the surface. 

Fig. 5 shows the dependence of the horizontal width of 
the precursor film w(t) on time for a typical simulation 
run for N = 1525. Qualitatively, the data look similar to 
that of Refs. (l4|jl(| with two regimes visible. We have 
analyzed the data as follows. First we assume that 

/ A(t-h) a +w , forti <t<t 2 , . . 
w{t) = \B(t-t 2 r+ W i, fort>t 2 . < 14 > 

in which t\ and t 2 denote cross-over times, wq and w' 
are constants, and a is the exponent for the first and 
for the second regime. ^From the data one expects to 
find two different power law regimes. A standard trick 
is to estimate the initial transient time t\ and then plot 
In («;(£) — wq) vs. In t. When trying to obtain (3 in a 
similar manner, the difficulty lies in the fact that the 
estimation of t 2 which denotes crossover towards final 
2D diffusion, is not as straightforward as for t\. A linear 
least-squares fit for the first regime gives w(t) ~ ^o.8±0.i 
The same procedure is then applied to the other regime. 



The slope of the least-squares linear fit for the second 
"diffusive" regime gives w(t) ~ t°- 5±01 . 

These two exponents obtained from the data are con- 
sistent with the results of Refs. HJlfJ. They found a 
crossover for the width of the precursor film w(t) from 
"almost linear" (~ t 9 ) to "diffusive" (~ t ' 5 ) behavior. 
When we translate this to the 3D situation, we recover 
the two "diffusive" regions with different effective trans- 
port coefficients in accord with experiments ]l[-p|,p^[ . Us- 
ing our result a ~ 0.8 and extracting A we find that the 
early-time diffusion coefficient D e « 5.4 x 10 _6 m 2 /s. 
For the late time "diffusion" coefficient we find Dg » 
1.2 x 10~ 6 m 2 /s and thus D e /Di sa 5. These values are 
somewhat larger than the measured ones that range be- 
tween 0.4- 2.0 x 1CT 12 m 2 /s f|. Also, typical experimen- 
tal ratios are of the order of 100 — 1000. The difference is 
not surprising since in our units T « 1600 K, and there 
are no surface diffusion barriers. Corrugation of the sur- 
face would most likely tend to lower the late time "dif- 
fusion" coefficient thereby making the ratio larger. We 
note that extrapolating our values of D to room temper- 
atures gives about 10 -15 m 2 /s, in reasonable agreement 
with experiments. 

The fact that w is a function of both time t and the 
number of dimers N suggests that there mi ght exist a 
scaling form for w(t, N) as suggested in Ref. jig] , which 
is of the following form: 

w(t) = t x <S>(t/N y ) . (15) 

However, for the present case we find no such scaling for 
the range of times and system sizes studied. This is in 
part because for our relatively small systems, crossover 
to diffusive behavior is very sharp and thus the data do 
not collapse. 

Shifted case 

For the shifted case, Figs. 6(a) and (b) show a typical 
evolution of the droplet for N = 1525 shifted dimers dur- 
ing spreading. Initial configuration is identical to that of 
the ordinary case. Again, the dimers that have not yet 
come into contact with the surface behave in a manner 
similar to the ordinary case. The molecules in the middle 
of the droplet are mainly oriented perpendicular to the 
surface. A striking difference between the ordinary and 
the shifted case is the development of a compact precur- 
sor layer which appears very well ordered at all stages 
of spreading. Fig. 7 shows the pair correlation function 
within the precursor film for the shifted case taken at 
t = 80000 r.u.. Clear peaks can be observed correspond- 
ing up to about fourth or fifth nearest neighbor dimers. 
The precursor layer in this case indicates a high degree 
of local ordering even at these elevated temperatures. 

We have also calculated the distribution of orientations 
for neighbouring dimers for the shifted case. The overall 
shape of the profile, which is shown in Fig. 8 corre- 
sponding to late stages of spreading, reflects the proper- 
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ties discussed above. A clear peak can be seen indicating 
orientation of nearby dimers in a preferred (vertical) di- 
rection. In this case the orientational correlation of the 
dimers extends through the entire precursor layer. 

Fig. 9(a) shows typical density profiles of dimers taken 
at three different times. Whereas in the ordinary case 
a step developed at the edges of the droplet, in the 
shifted case the edge of the precursor film always re- 
mains very sharp and well-defined. It is interesting to 
compare the calculated density profiles with the experi- 
mental ones shown in Fig. 9(b) for triloxane polyoxyethy- 
lene molecules spreading on silica bearing a dense grafted 
layer of trimethyls ||. In this case grafting results in a 
hydrophobic surface. The triloxane polyoxyethylene is a 
hammer-shaped molecule which has a hydrophobic (the 
trisiloxane head) and a hydrophilic (polyoxyethylene tail) 
part. The attraction of the hydrophobic group to the sur- 
face and the repulsion between the hydrophilic part and 
the surface forces the molecules lie perpendicular to the 
surface. Despite the enormous difference in the horizon- 
tal scales, the simulated and experimental profiles are in 
good qualitative agreement. 

Fig. 10 shows the dependence of w(t) for the shifted 
case. Again a crossover form "almost linear" to "diffu- 
sive" behavior can be seen. The crossover appears to 
be somewhat sharper in this case. With the smallest 
system size studied the spreading stops completely due 
to finite-size effects. With larger systems the spreading 
continues in a diffusive manner. We have analyzed the 
data for w(t) in the same way as before, and find that 
within the "almost linear" regime w(t) ~ £ -8±0-! anc i m 
the "diffusive" regime w(t) ~ t aA±0A . We again convert 
these results into the 3D case and recover the two dif- 
fusive regions. We estimate that D e ks 5.4 x 10 _6 m 2 /s, 
and Dg w 1.1 x 10~ 7 m 2 /s which give D e /Dg « 50, in 
good agreement with experiments ||. 

The reason for De being about an order of magnitude 
smaller than in the ordinary case can be understood from 
simple energetic arguments. The activation energy for 
diffusion in the shifted case is larger due to the high de- 
gree or local ordering. We have estimated the activation 
energies E a for diffusion for the two cases by calculating 
the average energy due to neighbours of a dimer located 
within two bond lengths from the edges of the droplet. In 
the ordinary case we find E a w 0.9 eV and in the shifted 
case E' a w 1.2 eV. If we make the assumption that Dg fol- 
lows the Arrhenius form D( oc e ~ E "/ kT ; W e can estimate 
that 

e -EJkT 



This is fully consistent with results extracted from the 
width of the precursor film. 

In the case of shifted dimers we again checked the scal- 
ing form of Eq. (16), but did not find a good data collapse 
for the present times and system sizes studied. 



B. Tetramers 

For chains consisting of four monomers we have studied 
three different cases, namely two ordinary cases (rather 
stiff and completely flexible tetramers) and a completely 
flexible shifted case. The surface potential for the grafted 
end is set to be Vj and for the ordinary cases we employ 
the surface potential Vjy. In the shifted case the sur- 
face potential of the other end is set to Vy, whereas the 
monomers in the middle of the chain do not have any 
interactions with the surface. 

Rather stiff tetramers 

Rather stiff tetramers are characterised by a fairly 
strong angle-dependent potential between consecutive 
bonds in a chain, namely Vq = eg (cos # + 1), where 
eg = 10e/. Figs. 11(a) and (b) show the evolution 
of a droplet for N — 785 in a typical simulation run. 
The chains appear to be initially fairly straight and the 
calculation of the distribution bond angles confirms this 
observation. The apparent holes in the droplet are due 
to thermal fluctuations. The calculation of distribution 
of orientations for neighbouring chains also supports the 
conclusion that the structure of the initial configuration 
is disordered and liquid-like. 

Qualitatively, the results are similar to the case of ordi- 
nary dimers. The chains in the middle of the droplet are 
oriented approximately perpendicular to it, the grafted 
end being closest to the surface. Closer to the edges, the 
chains tend to attain a more horizontal orientation. The 
centre of the droplet appears to have a very complicated 
structure. This is due to the finite flexibility of chains, 
as calculation of the bond angle distribution reveals that 
the chains are bent slightly more than in the initial con- 
figuration. 

We have also calculated the pair correlation function 
within the precursor layer corresponding to late stages 
of spreading. The structure of the precursor film in this 
case is disordered and liquid-like. At very late stages the 
neighboring chains tend to orient themselves in the same 
direction. This is revealed by calculating the distribution 
of orientations. The final stages of spreading correspond 
to diffusively thinning monolayer where the chains prefer 
to be relatively straight. 

The density profiles of Fig. 12 are characterised by a 
rounded overall shape with a peak corresponding to the 
chains that are not yet in the precursor film. This is 
consistent with the observations made from the configu- 
rations of Fig. 11. Fig. 13 shows w(t) for three different 
system sizes. A noteworthy feature is that the crossover 
towards late time diffusive behavior is not as clear as in 
the case of dimers. For the largest system size we find 
that w(t) ~ t° 01 within the initial "almost linear" 
regime, but our data for this system size does not extend 
far enough to capture the second "diffusive" regime. For 
the smallest system size N — 505 we find that in the "al- 
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most linear" regime w(t) ~ £°- 8±al and in the "diffusive" 
regime w(t) ~ t°- 5±Q1 . ^From these data we estimate 
that D e w 5.0 x f0" 6 m 2 /s and w 7.0 x 10" 7 m 2 /s 
which give D e /De ~ 7. 

For the present case, we also find a good data collapse 
using the scaling form of Eq. (16), with x = 0.9 and 
y = 0.9. Fig. 14 shows the scaling function $(z) ~ const. 
for z < 1, and $(2) - z 1/2_I for z > 1 @. 

Completely flexible tetramers 

The results for completely flexible tetramers are very 
similar to the results for rather stiff ones. However, in 
this case the density profiles in Fig. 15 appear to be some- 
what flattened at late stages as compared to the rather 
stiff case. This is evidently due to the greater flexibility of 
the chains. The width of the precursor film for N = 505 
tetramers gives w(t) ~ ^0.9±o.i m ^he "almost linear" 
regime, and w(t) ~ ^o.5±o.i m ^ e "diff us i ve " regime. We 
find D e w 4.5 x 10" 6 m 2 /s and D e w 2.6 x 10" 7 m 2 /s 
which give D e /D( « 20. Scaling of Eq. (16) is again 
obeyed, with x — 0.9 and y — 0.9. 

Completely flexible tetramers shifted case 

For the shifted tetramer case, the grafted end has the 
usual surface potential Vj whereas the other end had 
potential Vy for which the equilibrium distance from 
the surface extends about three bond lengths further 
away. We thus expect the results to be similar to those 
for the shifted dimers, with possible differences arising 
from the greater configurational entropy of the chainlikc 
molecules. Here we mostly concentrate in the morphol- 
ogy of the droplets during spreading. 

Figs. 16(a) and (b) show the evolution of the droplet 
for N — 1010. It can be seen that the evolution is strik- 
ingly different from the ordinary case. As in the case 
of shifted dimers, the precursor films appears to be very 
compact and well-ordered. We have calculated the pair 
correlation function within the precursor film, which is 
shown in Fig. 17. Clear peaks corresponding up to 
about fourth or fifth nearest neighbour chains are clearly 
present, which is an indication of a high degree of local 
ordering. If one compares the pair correlation functions 
between shifted dimers and tetramers, it can be seen that 
for tetramers the peaks appear to be somewhat broad- 
ened. This is evidently due to the chain-like structure of 
the tetramers. Thus, with the present choice of interac- 
tions the influence of chain flexibility is rather small even 
at high temperatures. 

We have also followed the time evolution of the den- 
sity profile of the droplet. This is shown in Fig. 18. The 
profile develops from a fairly sharply peaked one towards 
a smooth but compact form. These profiles again con- 
firm the observations made from the snapshots. We have 
estimated the effective diffusion barriers at late times 
and find that E' a ~ 2.9 eV compared to w 1.2 eV for the 



shifted dimer case. This can be understood on the basis 
that within a well ordered layer, the number of neigh- 
bors for tetramers should be roughly more than twice 
the corresponding number for dimers. 

C. Octamers 

The results for chains built up of eight monomers (oc- 
tamers) are presented in this section. We have studied 
two different systems, namely one in which the chains 
are rather stiff and another in which the chains are com- 
pletely flexible. The grafted end has the usual surface 
potential Vj whereas the other monomers have Vjy. It 
should be pointed out that due to CPU-time constraints 
out system sizes are relatively small. 

Rather stiff octamers 

Initial configuration for a droplet of rather stiff oc- 
tamers (ee = 10e/) of size N = 488 is shown in Fig. 
19(a). It is characterised by a fairly complex structure. 
The tendency of the chains to be relatively straight is 
clearly visible. Fig. 19(b) shows a typical evolution of the 
droplet. The structure in the middle of the precursor film 
is very complex. Chains at the edges of the droplet again 
tend to lie parallel to the surface. We have calculated 
different time-dependent quantities such as the pair cor- 
relation function within the precursor film taken at late 
stages of spreading. From the shape of the function we 
can immediately conclude that the precursor film is dis- 
ordered and liquid-like. The directions of neighbouring 
chains are fairly strongly correlated, however. Calcula- 
tion of the distribution of bond angles reveals that the 
chains tend to remain fairly straight through the whole 
spreading process. 

A set of typical density profiles are shown in Fig. 
20. It can be seen that at later times they become 
rather smooth and rounded. For the width of the pre- 
cursor film we find that initially w(t) ~ f0.9±o.i with 
D e w 1.7 x 10~ 5 m 2 /s. For the late-time behavior, we 
find that w(t) ~ <°- 5 ± 01 and D e w 4.3 x 10" 6 m 2 /s. For 
the ratio we thus find D e j Di w 4. As far as the scaling 
is concerned, we were not able to collapse the data for 
different system sizes for the present case. 

Completely flexible octamers 

The results for N — 488 completely flexible octamers 
are very similar to results for the rather stiff ones. Again 
the precursor film is disordered and a crossover from 
"almost linear" (w(t) ~ ^o.8±o.i^ towards "diffusive" 
(w(t) ~ £0.5±o.i^ behavior for w(t) is recovered, with 
D e « 3.7 x 10" 5 m 2 /s and D e « 1.77 x 10" 6 m 2 /s, 
which yield D e /Di ps 20. The density profiles shown 
in Fig. 21 bear a close resemblance to the ones for rather 



7 



stiff octamers, the late-time profiles being somewhat less 
rounded in this case. 



D. The influence of Langevin dynamics 

We have also implemented Brownian dynamics into our 
computer code, the motivation being to study the influ- 
ence of a local thermostat for the present case. Recently, 
the NH thermostat has been claimed to be physically un- 
suitable for microscopic studies of droplet spreading [jl9| , 
despite the fact that for a small system out of equilib- 
rium, there is no unique "best" choice. We note that 
since we have a smooth surface, coupling to an auxiliary 
thermostat must be used here. 

To check our results, we have performed additional 
simulations for N = 555 shifted dimers 0. We have 
employed two different values for the friction coefficient 
?7, namely r]i = 3 x 10 14 s _1 and 772 = 0.3 x 10 14 s _1 . 
We have set our bare time step to 0.01 x 10~ 14 s. For 
771 we recover the two different regimes for w(t) with 
w(t) ~ i 10±01 and w(t) ~ t°- 5±0A . We can again extract 
the associated "diffusion" coefficients with the result that 
D e w 5 x 10" 6 m 2 /s and D e w 1.5 x 10~ 7 m 2 /s for the 
early and late time regime, respectively. For the ratio we 
thus find D e /Dg « 30. Calculation of the pair correla- 
tion function within the precursor film again reveals that 
the layer has a high degree of local order. 

For 772 we also recover the two regimes with w(t) ~ 
£0.9±o.i an( j w ^ ^ £0.4±o.i^ -p Qr £ n j g p ar ticular case 

D e w 8 x 10" 5 m 2 /s and D e « 4.0 x 10~ 7 m 2 /s for the 
early and late time regimes, respectively. For the ratio 
we obtain D e /Di « 200. Reducing the friction would 
further increase this ratio since the late-time regime is 
dominated by effective diffusion barriers that are inde- 
pendent of 77. 

Based on our additional results we can conclude that 
qualitatively and quantitatively our results and conclu- 
sions are unaffected by the choice of the thermostat. The 
main effect of the local thermostat with respect to the NS 
thermostat is a slightly smoother crossover towards the 
late-time regime. 

IV. CONCLUSIONS AND DISCUSSION 

In this work we have studied a simple model of dy- 
namics of spreading for rigid and flexible molecules that 
interact asymmetrically with respect to a solid surface. 
We have studied two different cases. In the ordinary case, 
the end potentials are of different strength, but the equi- 
librium position of the molecules on the surface is hori- 
zontal. In the shifted case, however, the other end of the 
molecule has an equilibrium distance that is compatible 
with the length of the chain, i.e. the equilibrium position 
is vertical with respect to the surface. The latter case in 
particular can be considered as an effective model for the 



case of hydrophobic and hydrophilic surface interactions 
1- 

One of our main results is that the microscopic struc- 
ture in the precursor film drastically depends on the na- 
ture of the asymmetrical interactions. For the shifted 
case, there is a high degree of local order present which 
makes the density profile of the droplet unusually sharp 
and flat. This result is in good qualitative agreement 
with a recent experiment on a physically similar sys- 
tem ||. Moreover, our model recovers the overall t 5 
time dependence of the radius of the precursor film, and 
we have been able to quantitatively estimate the associ- 
ated transport coefficients. Typical ratios of the early- 
time coefficients D e to the late-time ones Di are in very 
good agreement with the experiments. Furthermore, our 
model demonstrates how this ratio increases with increas- 
ing local order in the precursor film, in cases where the 
late-time diffusion is controlled by energy barriers arising 
from neighboring molecules. 

Recently, the choice of the global NH thermostat used 
here and in Refs. |l4|jl6| has been criticized by De Con- 
inck et al. |l9| ] on basis of the argument that in an inho- 
mogeneous system a global thermostat is not physically 
justified. However, due to our smooth surface the heat 
must dissipated by other means than coupling to sub- 
strate atoms held at constant temperature, as was done 
111 Ref. @. Moreover, it is a well known fact that there 
is no unique way of controlling the instantaneous temper- 
ature of a non-equilibrium system. We have performed 
our simulations mainly with the NH thermostat, but test 
runs with Brownian dynamics reveal that the same qual- 
itative and quantitative behavior persists. Moreover, the 
results of Ref. jl6| as well as those of the present work 
compare very well with the simulations of De Coninck 
et al. |IT|] , once differences in the geometries (cylindrical 
vs. spherical) are properly accounted for. De Coninck et 
al. found that for largest droplets the number of atoms 
in the first (i.e. the precursor) layer was well described 
by N(t) ~ 7j°- 85 ±0- 05 and the corresponding radius by 

R 2 (t) ~ /j0.82±0.06 g|] Q n the Qther handj in Qur ge _ 

ometry w(t) ~ N(t) ~ /j°- 9 ± 01 for the "almost linear" 
regime. The flux of particles into the precursor layer is 
the same in the two independent studies and therefore 
the results are equivalent. The fact that the slower late- 
time "diffusive" regime was not reported in Ref. [ p"9[ is 
probably due to insufficient simulation times, since their 
system sizes were rather large. Thus both the qualitative 
and quantitative features of the spreading phenomenon 
are fairly insensitive to the choice of thermostat as well 
as the geometry. 

To summarize, we hope to have further demonstrated 
in this work that the spreading phenomenon at micro- 
scopic length-scales is a very complicated process. It 
seems highly unlikely that the properties of all the dif- 
ferent cases studied here and in other works could be 
obtained from a more general framework. On the other 
hand, there are many features of spreading, such as the 
time-dependence of the precursor radius that are rather 
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insensitive to the details of interactions, or molecular 
structure of the liquid. This work calls for more system- 
atic and controlled experimental as well as theoretical 
work in order to further classify the properties of tiny 
liquid droplets spreading on a solid surface. 
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Surface interaction parameters i 




cn 


Symbol 


1.0 e f 


5.0 a/ 


Vl 


0.06 t s 


5.0 a, 


Vu 


0.02 e f 


7.3 a f 


Vm 


0.006 e/ 


5.0 a, 




0.01 e/ 


8.0 a s 





TABLE I. Surface interaction parameters and their sym- 
bols used in this study. 



FIGURE CAPTIONS 



Fig. 1. Schematic illustration of the different surface 
potentials used in this study. The "grafted" end always 
has potential V\. 



Fig. 2. (a) Initial configuration for the ordinary case 
of N = 1525 dimers. The grafted end is represented 
by a large filled circle, (b) Same system taken at t = 
30000 r.u.. 



Fig. 3. Pair correlation function within the precursor 
film at t = 80000 r.u.. Note the disordered and liquid-like 
structure of the film. 



Fig. 4. Density profiles for the case of ordinary dimers. 
The late-time shoulders are due to the dimers that lie flat 
on the surface. These and all the other density profiles 
have been smoothed to remove noise. 



Fig. 5. The width of the precursor film w(t) for the 
ordinary dimer case. It can be characterised by an "al- 
most linear" regime which crosses over to a "diffusive" 
one. 



Fig. 6. (a) Initial configuration for the shifted case 
of N = 1525 dimers. (b) Same system taken at t — 
30000 r.u.. 



Fig. 7. Pair correlation function within the precur- 
sor film at t = 80000 r.u.. Clear peaks can be observed 
corresponding up to fourth or fifth nearest neighbour, 
indicating that the film displays a high degree of local 
order. 



Fig. 8. Distribution of orientations d Q for the shifted 
dimer case taken at t — 80000 r.u.. Orientations of dimers 
are correlated over the whole system. 



Fig. 9. (a) Density profiles for the shifted dimer 
case, (b) Experimental density profiles for triloxane 



polyoxyethylene molecules (hydrophobic head with a hy- 
drophilic tail) spreading on silica bearing a dense grafted 
layer of trimethyls M. The profiles are strikingly similar. 



Fig. 10. w(t) for the shifted dimer case. It can be 
characterised by an "almost linear" regime which crosses 
over to "diffusive" one. The late-time "diffusion" coeffi- 
cient is roughly one order of magnitude smaller than for 
the ordinary case. 



Fig. 11. (a) The initial configuration for N = 785 
rather stiff tetramers. (b) Same system taken at t = 
80000 r.u.. Notice the complicated and disordered struc- 
ture of the precursor filmi due to the flexibility of the 
chains. 



Fig. 12. Density profiles for rather stiff tetramers 
taken at three different times. The profiles are fairly 
smooth and rounded. 



Fig. 13. w(t) for rather stiff tetramers with three 
different system sizes. 



Fig. 14. Scaled data for w(t) for the three different 
system sizes with x = 0.9 and y = 0.9. 



Fig. 15. Density profiles for completely flexible 
tetramers taken at three different times. The profiles 
are somewhat flatter as compared to the rather stiff case 
due to the flexibility of the chains. 



Fig. 16. (a) Initial configuration for the shifted case 
of N — 1010 completely flexible tetramers. (b) Same 
system taken at t = 80000 r.u.. Notice the appearance of 
a compact precursor film with sharp edges. 



Fig. 17. Pair correlation function within the precur- 
sor film at t — 80000 r.u.. Clear peaks can be observed 
corresponding up to fourth or fifth nearest neighbour, 
indicating that the film displays a high degree of local 
order. 



Fig. 18. Density profiles for the shifted case of com- 
pletely flexible tetramers taken at three different times. 
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Fig. 19. (a) Initial configuration for the case of N = 
488 rather stiff octamers. (b) Same system taken at t = 
80000 r.u.. The structure of the precursor film is again 
disordered and liquid-like. 



Fig. 20. Density profiles for the case of rather stiff 
octamers taken at three different times. The profiles are 
fairly smooth and rounded. 



Fig. 21. Density profiles for the case of completely 
flexible octamers taken at three different times. The pro- 
files are somewhat flatter than in the case of rather stiff 
octamers. 



12 



